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The dynamics of precessing binary black holes (BBHs) in the post-Newtonian regime has a strong 
timescale hierarchy: the orbital timescale is very short compared to the spin-precession timescale 
which, in turn, is much shorter than the radiation-reaction timescale on which the orbit is shrinking 
due to gravitational-wave emission. We exploit this timescale hierarchy to develop a multi-scale 
analysis of BBH dynamics elaborating on the analysis of Kesden et al. [1]. We solve the spin- 
precession equations analytically on the precession time and then implement a quasi-adiabatic ap¬ 
proach to evolve these solutions on the longer radiation-reaction time. This procedure leads to an 
innovative “precession-averaged” post-Newtonian approach to studying precessing BBHs. We use 
our new solutions to classify BBH spin precession into three distinct morphologies, then investigate 
phase transitions between these morphologies as BBHs inspiral. These precession-averaged post- 
Newtonian inspirals can be efficiently calculated from arbitrarily large separations, thus making 
progress towards bridging the gap between astrophysics and numerical relativity. 

PACS numbers: 04.25.dg, 04.70.Bw, 04.30.-w 


I. INTRODUCTION 

Observations suggest that astrophysical black holes are 
generally spinning mu and can form binary systems 
[5] . Spinning binary black holes (BBHs) are a promising 
source of gravitational waves (GWs) [SHS] for current and 
future detectors [SHIS]. BBH dynamics is remarkably 
complex and interesting, especially when both BBHs are 
spinning. BBH systems have three angular momenta, the 
two spins and the orbital angular momentum, all coupled 
to each other. Spin-orbit and spin-spin couplings cause 
these angular momenta to process, changing their orien¬ 
tation in space on the precession timescale [mils]. On 
the longer radiation-reaction timescale, GWs take energy 
and momentum out of the system, thus shrinking the or¬ 
bit [T9ll2Q]. These emitted GWs encode all the richness of 
the processional dynamics but are also more challenging 
to detect and characterize than GWs emitted by non- 
precessing systems pikSS] . 

Expanding on the analysis in our previous paper [T], 
we introduce a multi-timescale analysis of the dynamics 
of spinning, precessing BBHs. Multi-timescale analyses 
are commonly used in binary dynamics. For example. 


* d.gerosa@damtp.cam.ac.uk 
t kesden@utdallas.edu 
t u.sperliake@damtp.cam.ac.uk 
§ eberti@olemiss.edu 
^ rossma@rit.edu 


in eccentric binaries the orbital period, periastron pre¬ 
cession, and radiation-reaction timescales usually differ 
by orders of magnitude; the dynamics of these systems 
can be studied using techniques that explicitly exploit 
this timescale hierarchy, such as osculating orbital ele¬ 
ments [29] or the variation of constants [30]. Exploit¬ 
ing timescale hierarchies leads to deeper understanding 
of the dynamics because different physical processes are 
decoupled and individually analyzed. 

Precessing BBHs evolve on three distinct timescales: 

1. BBHs orbit each other (changing the binary sepa¬ 
ration r) on the orbital time torb ~ ^^/^/(GM)^/^, 

2. the spins and the orbital angular momen¬ 
tum change direction on the precession time 
Ve - 

3. the magnitudes of the orbital energy and angular 
momentum decrease on the radiation-reaction time 
tRR c^rV{GM f. 

Here r = |r| is the magnitude of the binary separation, 
M is the total mass of the binary, and prefactors of order 
unity have been omitted. In the post-Newtonian (PN) 
regime, r » GMjc? and these timescales are widely sep¬ 
arated: 

forb ^ ^pre < fRR ■ (1) 

BBHs complete many orbits before their angular mo¬ 
menta appreciably precess, and the angular momenta 
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complete many precession cycles before the separation 
decreases significantly. 

The first inequality (torb ‘C tpre) has been widely ex¬ 
ploited to understand spin dynamics and approximate 
the GW signal. This approximation forms the founda¬ 
tion of the orbit-averaged spin-precession equations for 
adiabatic quasicircular orbits examined extensively in the 
pioneering study of Apostolatos et al. m and later ex¬ 
tended by Arun et al. |3ll|32]. Using these equations, 
several authors have systematically explored the physics 
of spin precession and their implications for GW detec¬ 
tion [SSI [S3] and astrophysics [SS1|SS]. Following the early 
work by Schnittman on spin-orbit resonances [37], PN 
spin dynamics has been used to predict the final spins [SS] 
and recoils [35] [SB] following BBH mergers, select initial 
conditions for numerical-relativity simulations |39j , char¬ 
acterize formation scenarios for stellar-mass BH binaries 
m, and address the distinguishability of these scenarios 
by future GW observations [STHiS] . 

The second inequality (tp^e ^ t^n) has received less 
attention because until now there were no explicit so¬ 
lutions for generic spin precession (unlike the Keplerian 
orbits that readily allowed orbit averaging in previous 
work). Our new solutions for spin precession allow us to 
fully exploit the timescale hierarchy of Eq. Q , expanding 
and detailing the ideas put forward in our previous Let¬ 
ter [1]. We showed that spin precession is quasi-periodic 
implying that the relative orientations of the three an¬ 
gular momenta are fully specified by a single parame¬ 
ter, the magnitude of the total spin, that oscillates on 
the precession time. As is common in multi-timescale 
analyses, once the dynamics on the shorter time has 
been solved, the behavior of the system on the longer 
time scale can be studied as a quasi-adiabatic process. 
We evolve our precessional solutions during the inspi¬ 
ral by double-averaging the PN equations over both the 
orbital and the precessional timescales. Semi-analytical 
precession-averaged inspirals turn out to be extremely 
efficient and can be carried out from infinitely large sep¬ 
aration with negligible computational cost. 

While our focus in this work is on spin precession, our 
study benefits from several recent investigations which 
also used separation of timescales to efficiently and accu¬ 
rately approximate both the dynamics and the associated 
GW signal. A series of papers by Klein et al. [zasniis] 
used a multi-scale analysis to construct semianalytic 
approximations to the frequency-domain waveforms for 
generic two-spin processing binaries. Lundgren and 
O’Shaughnessy [4^ used this timescale hierarchy to con¬ 
struct semianalytic approximations to the inspiral of pro¬ 
cessing binaries with a single significant spin. The GWs 
emitted during the full inspiral-merger-ringdown of spin¬ 
ning, processing binaries were also investigated using phe¬ 
nomenological models based on a single “effective spin” 
approximation [STMS] and the effective-one-body frame¬ 
work [50] . 

This paper is organized as follows. In Section [IT] we de¬ 
rive explicit solutions for generic BBH spin precession at 


2PN order on timescales short compared to the radiation- 
reaction time tRR. These solutions allow spin precession 
to be classified into three different morphologies charac¬ 
terized by the qualitative behavior of the angle between 
the components of the two spins in the orbital plane. In 
Sec. |III[ we use our new solutions to precession average 
the radiation reaction on the binary at IPN order and 
demonstrate how this precession averaging improves the 
computational efficiency with which GW-induced inspi¬ 
rals can be calculated compared to previous approaches 
relying solely on orbit averaging. Precession-averaged 
evolution does not preserve the memory of the initial pre¬ 
cessional phase, just like orbit-averaged PN evolutions 
do not track the orbital phase. Sec. [IV] explores phase 
transitions between the three precessional morphologies, 
which are readily identified using our new formalism and 
have potentially interesting observational consequences. 
Finally, we conclude in Sec.jVj highlighting the relevance 
of our new PN approach to both the theoretical under¬ 
standing of BBHs and observational GW astronomy. We 
mainly focus on the relative orientation of the momenta; 
the evolution of the global orientation of the system will 
be addressed elsewhere m- Throughout the paper, we 
use geometrical units (G = c = I) and write vectors in 
boldface, denoting the corresponding unit vectors by hats 
and their magnitude as (e.g.) L = |L|. Latin subscripts 
(j = 1, 2) label the BHs in the binary. Binaries are stud¬ 
ied at separations r > lOM, taken as a simple but ad hoc 
threshold for the breakdown of the PN approximation 
[5SII54j . Animated versions of some figures are available 
online at the URLs listed in Ref. [55] . 


II. ANALYTIC SOLUTIONS ON THE 
PRECESSIONAL TIME SCALE 


In this section we focus on the binary dynamics on 
the precessional time. Angular momentum conservation 
(Sec, jll A ) and the existence of a further constant of mo¬ 
tion (Sec. IIBI provide a simple parametrization of the 
binary dynamics through the identification of effective 
potentials. Solutions are classified according to the pre¬ 
cession geometry (Sec. 


an inertial frame (Sec. 


IIDl. 


IIG) and eventually expressed in 


A. Parametrization of precessional dynamics 

Let us consider BBHs on a circular orbit Q Let mi 
and m 2 denote the BBH masses, in terms of which we 


^ Our approach can be readily generalized to nonzero eccentricity 
without complicating the geometry since the orbital pericenter 
processes on a shorter timescale than the BBH spins do. We 
restrict our attention to circular orbits since radiation reaction 
is expected to suppress the eccentricity at large separations for 
most astrophysical systems [I9l[20]. 
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FIG. 1. Reference frames used in this paper to study BBH 
spin precession. The angles 6 \, 62 , A<1>, and 612 are defined in 
a frame aligned with the orbital angular momentum L (left 
panel). The binary dynamics can also be studied in a frame 
aligned with the total angular momentum J (right panel). 
Once L is taken to lie in the a; 2 -plane, its direction is spec¬ 
ified by S through the angle 9l- The angle ip' corresponds 
to rotations of Si and S 2 about the total spin S. The two 
frames pictured here are not inertial because the direction of 
L changes together with the spins to conserve J. These angles 
are defined in Eqs. and §. 


can define the mass ratio q = milTni < 1, the to¬ 
tal mass M = mi -|- m 2 , and the symmetric mass ra¬ 
tio r] = mim 2 /M^. The spin magnitudes Si = mfxi 
(i = 1,2) are most conveniently parametrized in terms 
of the dimensionless Kerr parameter 0 < Xi < 1, while 
the magnitude of the orbital angular momentum is re¬ 
lated to the binary separation r through the Newtonian 
expression L = ri{rM^Y^^. 

The three angular momenta L, Si and S 2 in principle 
constitute a nine-dimensional parameter space. However, 
there exist numerous constraints on the evolution of these 
parameters, greatly reducing the number of degrees of 
freedom. At the PN order considered here, the magni¬ 
tudes of both spins are conserved throughout the inspi¬ 
ral (see e.g. Ref. |12]), reducing the number of degrees 
of freedom from nine to seven. The magnitude of the 
orbital angular momentum is conserved on the preces¬ 
sion time (although it shrinks on the radiation-reaction 
time), further reducing the number of degrees of free¬ 
dom from seven to six. The total angular momentum 
J = L-I-S 1 -I-S 2 is also conserved on the precession time, 
reducing the number of degrees of freedom from six to 
three. As described in greater detail in the next subsec¬ 
tion, the projected effective spin ^ [56l El] is also con¬ 
served by both the orbit-averaged spin-precession equa¬ 
tions at 2PN and radiation reaction at 2.5 PN order, pro¬ 
viding a final constraint that reduces the system to just 
two degrees of freedom. In an appropriately chosen non- 
inertial reference frame precessing about J, precessional 
motion associated with one of these degrees of freedom 


can be suppressed, implying that the relative orienta¬ 
tions of the three angular momenta L, Si and S 2 can 
be specified by just a single coordinate! We will provide 
an explicit analytic construction of this procedure in this 
and the following subsection. 

We begin by introducing two alternative reference 
frames in which the relative orientations of the three an¬ 
gular momenta can be specified explicitly. As shown in 
the left panel of Fig. one may choose the z'-axis to lie 
along L, the cc'-axis such that Si lies in the a;'z'-plane, 
and the y'-axis to complete the orthonormal triad. In 
this frame only three independent coordinates are needed 
to describe the relative orientations of the angular mo- 


menta; we choose them to be the angles 


cos 01 = Si • L , 


(2a) 

cos 02 = S 2 • L, 


(2b) 

. Si X L 

S 2 X L 


COSA$ = ' 


(2c) 

|Si xLj 

IS 2 X L|’ 



where the sign of A$ is given by (cf. Fig.[^ 

sgn A$ = sgn{L • [(Si x L) x (S 2 x L)]}. (2d) 


The relative orientations of the three angular momenta 
can alternatively be specified in a frame aligned with the 
total angular momentum J. For fixed values of L, Si, 
and S 2 , the allowed range for J = jJj is 


Jmin Si J Si Jn 


(3a) 


where 


Jmin = max(0, L - - S 2 , [S'! - S '21 - L) , (3b) 

Anax = L + Si + S 2 ■ (3c) 

As shown in the right panel of Fig. one can choose the 
z-axis parallel to J and the cc-axis such that L lies in the 
ccz-plane: 

J = Jz and L = Lsin6*LX-b Tcos6*lz . (4) 

The third unit vector y = z x x completes the orthonor¬ 
mal triad. The total spin S = Si-|-S 2 = J — L will also 
lie in the xz-plane: 


S = —Lsind^x + {J — Lcos9l)z , 


implying 


cosOl = 


J2 -b - S'2 


2JL 

We can also define a unit vector 

~ {J — Lcos 6 l)x + LsiuOlz 

3 I = -- - -- 


(5) 

( 6 ) 

(7) 


which also lies in the xz-plane but is orthogonal to S. 

While the magnitudes L and J of the orbital and to¬ 
tal angular momenta are conserved on the precession 
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timescale, the same is not true for the total-spin mag¬ 
nitude S', which oscillates within the range 


<S <S„ 


-^min _ ^ _ 


(8a) 


where 


5’min = max(| J - L|, |5i - 52|) , (8b) 

>8*max — min( J-I-L, Si-I-S2) . (8c) 


S can be used as a generalized coordinate to specify 
the directions of the angular momenta J, L, and S; we 
can see from Eqs. Q-® that it is the only coordinate 
needed to specify these directions in the xyz-iiame. 

Specifying the directions of the individual spins Si and 
S2 in the xyz-irame requires an additional generalized 
coordinate, which can be chosen to be the angle if' be¬ 
tween Sx in Eq. (j^ and the projection of Si into the 
plane spanned by Sx and y, as shown in the right panel 
of Fig. This angle corresponds to rotations of Si and 
S2 about S and is given analytically by 


cos ip' 


Si-Sx 
|Si X S| ■ 


(9) 


In terms of the two coordinates S and ip' varying on the 
precession timescale, the three angular momenta in the 
xyz-iiame are 


^ A 1 A 2 . J^ + L^-S- 

L = _ ^ X -|- ---Z 


Si = 


2J 

52 + 


2J 


q2 q2 
Oi 


25' 


A 3 A 4 

25 


(10a) 

(cos yi'Sx -I- sin p'y) 


1 


r[-(52 + 5?-5|)^iA2 


4J52' 

-I- (+ S'^)A 3 A 4 cos 
+ ^^3^4 smtp y 


1 


S2 = 


25 


[(52 -h 51^ - 5|)(J2 - + 5-2) 

(10b) 

(cos yi'Sx -I- sin p'y) 


4J52 

-I- Ai A2A3A4 cos p']z , 


52-5?x >13414 


S - 


25 


4J52 


[(52 + 5| - Sf)A4A2 


+ (+ S‘^)A 3 A 4 cos p']k 
- ^A3^4sm(p y 


1 


[(5^ + Sl- Sl){J^ -L^ + 5^) 


4J52 

— Ai A2A3A4 cos p'\z , 
where we defined: 


(10c) 


= [j2 _(i_ 5)2]i/2^ (11a) 

>l2= [(L + 5)2 -j 2]1/2^ (lib) 

^3 = [5'-(5i- 52)2] 1/2, (11c) 

414= [(5i+52)2-52]1/2. (lid) 


All the Ai’s are real and positive in the ranges specified 
by Eqs. (|^ and 


B. Effective potentials and resonances 


As anticipated in the previous subsection, there is an 
additional conserved quantity that can be used to elimi¬ 
nate p' and thereby specify L, Si, and S2 with the single 
generalized coordinate 5. This quantity is the projected 
effective spin 


^ = M-2[(l + (7)Si + (l + g-i)S2]-L (12) 


which is a constant of motion of the orbit-averaged spin- 
precession equations at 2PN order and is also conserved 
by radiation reaction at 2.5 PN order. Using Eqs. (101, 
we can express ^ as a function of 5 and p' 


as, p') = {(J2 - l 2 _ sais^l + qf - (52 - 52)(1 - 92)] 

— (1 — 92)21^^2^3^4 cos:p'}/( 49 M^ 52 L). 

(13) 


Conservation of ^ restricts binary evolution to one¬ 
dimensional curves ^(5, tp') = ^ in the 5(/?'-plane as 
shown in the right panel of Fig. The simple depen¬ 
dence of as, 'f') on p' motivates us to define two “effec¬ 
tive potentials” [T] corresponding to the extreme cases 
cos:^' = =f 1 for which L, Si and S2 are all coplanar: 


e±(^) = -L^- sais^a+qf - (52 - 522)(i - 9")] 

± (1 - 9')AiA2A3A4}/(49M252l) . (14) 

At 5niin and 5niax 

^-{S^in) = C+(5mi„) , e-(5„,ax) = ^+{Sm.a , (15) 


because one of the A^’s defined in Eqs. (11) vanishes if 
5 = 5niin or 5 = 5max- The functions ^±XS) thus form 
a loop that encloses all allowed values of 5 and a as 
shown in the left panel of Fig. BBHs are constrained 
to evolve back and forth along horizontal line segments of 
constant ^ bounded by the two effective potentials ^±{S). 
The turning points in the evolution of 5 are given by 
the solutions of ^±{S) = a where the binary meets an 
effective potential. Once squared, the equation (5) = 
^ reduces to the following cubic equation in 5^: 


aeS^ + (745^ + CT252 + (To = 0, 


(16a) 
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FIG. 2. Left: Effective potentials C±(<S') for BBHs with q — 0.6, Xi = Xa = Ij ^ = lOOM, and J = 2.34M^. Conservation of the 
projected effective spin ^ constrains the BBH spins to process along horizontal lines bounded by the effective-potential curves. 
The horizontal dashed lines intersecting the effective potentials at Smin and Smax (marked by empty square^ divide BBH spin 
precession into three different morphological phases distinguished by whether the angle tp' defined by Eq. ® oscillates about 
TT (top orange region), circulates from 0 to 27r (middle grey region), or oscillates about 0 (bottom purple region). The effective 
potentials admit two extrema ^min and ^max (marked by empty tr iang les) corresponding to the spin-orbit resonances discovered 
in Ref. [^. Right: Contours of constant given by Eq. (131 for the same binary parameters. As BBH spins precess 

along the horizontal dashed lines in the left panel, they move along the curves in the ^(/p'-plane in the right panel illustrating 
the three morphological phases. 


where 

aG = q{l + q)‘^, (16b) 

(74 = (1 -|- q)'^[ — 2J^q ^1 -|- -|- 2LM^ 

-b (1 - g) (S'! - qSf )], (16c) 

a 2 = 2il + q)\l-q)[J^iqSf-Sl) 

-L^{S!-qSi)] + q{l + qr{J^-Ly 
- 2LM^f,q{l -b g)[(l -b q){J^ - L^) 

+ (1 - q){Sl - ^1)] + AL^M^eq" , (16d) 

ao = {l-q^)[L^{l-q^)iSf-Sir 
-(l + q)iqSf-Si)iJ^-Ly 
+ 2LM^qas!-Sl)iJ^-L^)], (16e) 

which admits at most three real solutions for 5* > 0. The 
number of solutions in the range allowed by Eqs. ([^ must 
be even because the two effective potentials form a closed 
loop and the Jordan curve theorem requires the number 
of intersections between a continuous closed loop and a 
line to be even (although these intersections can co¬ 
incide at extrema). Since two is the largest even number 
less than three, the equation ^±(5') = ^ will generally 
have two solutions which we denote by S± {S- < S+). 

The total-spin magnitude S will oscillate between S- 
and 5'+ implying that spin precession is regular or quasi- 


periodic (this will be shown explicitly in Sec. IID below). 
The motion of the spins is not fully periodic because in 
an inertial frame the basis vectors x and y will generally 
not precess about J by a rational multiple of tt radians 
in the time it takes S to complete a full cycle from S- 
and 5+ and back again. The turning points S = S± lie 
on the effective potentials, implying from the definition 
cos(p' = =f 1 that all three vectors L, Si, and S2 are 
coplanar. The qualitative evolution of ip' is related to 
the nature of the turning points S±. This is illustrated 
in Fig. 1^ where horizontal lines in the effective-potential 
diagram (left panel) correspond to contours of constant 
^{S,(p'), computed using Eq. (13) (right panel). Three 
different cases are possible. 


1. Both turning points lie on 

^+(5+)=e+(5-)=e- (17a) 

ip' oscillates about tt never reaching 0 (orange re¬ 
gion in Fig. [^. 

2. One turning point is on and the other is on 

C±(5-)=CT(^+) = ^ (17b) 

p' monotonically circulates from —tt to tt during 
each precession cycle (grey region in Fig.[^. 
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3. Both turning points lie on 


m 


e_(5+)=e-(5-)=^ (17c) 

ip' oscillates about 0 never reaching tt (purple region 
in Fig. [^. 


The boundaries between the three regions are given by 
those values of ^ at which one of the turning points S± co¬ 
incides with either S'min or iFmax (dashed lines in Fig.[^. 
Note that ^(S'min) may be less or greater than ^(S'max) 
depending on the values of g, Xi, r and J. 

The two turning points are degenerate {S+ = S'-) at 
the extrema ^min and ^max of the effective potentials. At 
these extrema the derivatives 


dS 


± 


1 + g 
2 qM‘^S^L 
1-9 


(1 _ g)( J2 _ - Si) - (1 + q)S^ 

A A A A s^ - + Sf + Sl)S^ 

A1A2A3A4L 12; 

-f (J2 + L^){Sl - SD'^S'^ + (Si + SDiJ^ - L'^)'^S^ 


- (sf - - Ly 


(18) 


vanish and S = S_ = S+ is constant. Since 


lim > lim 

S —>-S’min dS 5—VSmi 

lim < lim 

S — >Sraa.x dS S —^Sma 


dy 

dS ’ 
dS ’ 


(19a) 

(19b) 


and at most two turning points can exist, it follows that 
admits a single maximum in [Smin,Smax] and ad¬ 
mits a single minimum in [Smin, Smax]- The effective po¬ 
tentials therefore have exactly two distinct extrema for 
each value of the constants J, r, g, xi and X 2 - As clar¬ 
ified below, these special configurations correspond to 
the spin-orbit resonances discovered by other means in 
Ref. [37]. 

The equal-mass limit 9 —1 corresponds to ^+{S) = 
yS) [cf. Eq. @] implying that S is constant for all 
values of ^ [note that ^±(S'min) ^±(5'inax)]- This fact 

was noted at least as early as 2008 by Racine m and it 
was recently exploited in numerical-relativity simulations 
[39l [59] , but the constancy of S' is a peculiarity of the 
equal-mass case and does not hold for generic binaries. 


C. Morphological classification 


Although the evolution of ip' already provides a way to 
characterize the precessional dynamics (Fig. [2), a more 
intuitive understanding can be gained by switching back 
to the L-aligned frame illustrated in the left panel of 
Fig. Substituting Eqs. (10) and (13) into Eq. ([^, we 
can express the angles 0i, 82 and Ad) as functions of S, 
J and This yields the remarkably simple expressions 


cos 6*1 = 


cos 6*2 = 


cos Ad* = 


1 

2(1-9)^! . 
9 

2(1 - q)S 2 


J2 - 2 qM^i 

L 1 -|- q 


-^-^ 

L 1-tq 


cos 012 — COS 01 COS 02 
sin 01 sin 02 


( 20 a) 

( 20 b) 

( 20 c) 


where the angle 0 i 2 = arccosSi • S 2 between the two 
spins can also be written in terms of S: 


cos 012 = 


52 - Sf - Si 
2 S 1 S 2 


( 20 d) 


Equations (20) parametrize double-spin binary preces¬ 


sion using a single parameter S. Some examples of the 
evolution of these angles over a precessional cycle are 
given in Fig. [^ The evolution of 0i and 02 is monotonic 
as S evolves between its two turning points S±; over a 
full precessional cycle these angles oscillate between two 
extrema lying on the effective potentials (dotted curves 
in Fig. [^. The evolution of Ad) can be classified into 
three morphological phases similar to that of ip': 

1 . Ad* oscillates about 0 (never reaching tt) if 


Ad>(S'_) = Ad>(5+) =0; 


( 21 a) 


2. Ad* circulates through the full range [— 7 r, 7 r] if 

Ad)(S'±) = 0 and Ad)(S' 4 ;) = tt ; ( 21 b) 

3. Ad* oscillates about tt (never reaching 0) if 

Ad>(5_) = Ad>(S'+) = 7 r. (21c) 


The evolution of Ad? allows us to unambiguously cate¬ 
gorize the precessional dynamics into the three different 
classes listed above. We refer to these classes as mor¬ 
phologies because of the different shapes traced out by 
the BBH spins over a precession cycle. We show some 
examples of how the allowed region inside the effective- 
potential loop is divided between these three morpholo¬ 
gies in Fig. [^ BBHs in the two oscillating morphologies 
are adjacent to the extrema of the effective potentials 
(^min and Cmax), while circulating binaries (if present) fill 
the gap in between. Schnittman’s spin-orbit resonances 
m can be reinterpreted as the limits of the two oscil¬ 
lating morphologies when the “precessional amplitude” 
S'+ — S- goes to zero at ^min and ^max) much like how 
circular orbits are the limits of eccentric orbits as the 
amplitude of the radial oscillations goes to zero. 

According to the criteria listed in Eqs. (21), bound¬ 
aries between the three morphologies (shown by horizon¬ 
tal dashed lines in Fig. ® occur at values of ^ where 
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S/M'^ 

FIG. 3. Analytical solutions given by Eq. ( |20[ ) for the evo¬ 
lution of the angles 6 \ (top panel), 62 (middle panel), and 
A® (bottom panel) during a precession cycle. The evolution 
of three binaries with ^ = 0.25 (blue), 0.3 (green) and 0.35 
(red) is shown for q = 0.8, xi = 1) X 2 = 0.8, r = 20M and 
J = 1.29M^. The evolution of di and 62 is monotonic during 
each half of a precession cycle and is bounded by the dotted 
lines for which cosifi = [these curves can be found by sub¬ 
stituting C±('S') for ^ in Eq. ( |20[ )]. Three classes of solutions 
are possible and define the binary morphology: Ail? can oscil¬ 
late about 0 (^ = 0.25), circulate (^ = 0.3) or oscillate about 
IT — 0.35). An animated version of this figure is available 
online at Ref. [SS], where precession solutions are evolved on 
iRR. 


COS A$ given by Eq. (20c I changes discontinuously at one 


of the turning points S± along the effective-potential loop 
C± (S). We know that A<1> is either 0 or tt along (S') be¬ 
cause L, Si, and S2 are coplanar when cos:^' = ±1 (see 
Fig. 0. A discontinuity can only occur when the denom¬ 
inator of Eq. ( |20c[ ) vanishes, i.e. where one of the spins 
is either aligned or anti-aligned with the orbital angular 
momentum (sin^i = 0). These discontinuities can only 
happen at the turning points S± because of the mono¬ 
tonic evolution of 9i during each half of the precession 
cycle, as shown in the top and middle panels of Fig. 
The four contours cos 9i = ±1 (sin^i = 0) are shown by 
dotted curves in Fig. we see that a boundary between 
morphologies occurs whenever these curves are tangent to 
the effective-potential loop ^±{S). These boundaries had 
previously been described as unstable resonances E7] . 

The geometrical constraints imposed by Eqs. (1^ and 
([^ imply that some morphologies may not be allowed for 
given values of L,J,q,xi, and X 2 - Three qualitatively 
different scenarios can occur, exemplified by the three 
panels of Fig. 

1. Left panel: BBH spins process in all three of the 
morphologies listed in Eq. (21). Libration about 


the coplanar configuration Ad? = 0 occurs for val¬ 
ues of ^ close to libration about the Ad? = tt 
configuration is found near ^max, and Ad? circulates 
for intermediate values of Our analysis in Ref. [T] 
was restricted to this case. 


2. Middle panel: Ad? oscillates about tt for ^ close to 
both and ^max, with circulation still allowed 
for intermediate values of f. 

3. Right panel: Ad? oscillates about tt for all values 
■fniin < f < 'fmax (circulation and oscillation about 
0 are both forbidden). 

To distinguish these scenarios, it is useful to examine 
the values of Ad? on the effective-potential loop at the 
extrema ^min and ^max- Although it is straightforward 
to evaluate Ad? numerically at ^max) one can gain more 
intuition by instead evaluating it at S'min- The value of 
Ad? is the same at these two points since the slope of the 
effective-potential loop ^+{S) connecting them is positive 
while that of the cos^i = ±1 contours is negative (as can 
be seen in Fig.|^. The curves therefore cannot be tangent 
to each other implying that Ad? must remain constant on 
this portion of the loop. Equation ( |8b| ) requires that S'min 
equals the greater of \ J — L\ and |Si — S 2 I; in the former 
case L and J are anti-aligned, while in the latter case Si 
and S 2 are anti-aligned. In either case, the components of 
Si and S 2 perpendicular to L are anti-aligned (Ad? = tt). 
This implies that Ad? will oscillate about tt near ^max for 
all values of J, L, Si, and S 2 (as can be seen in all three 
panels of Fig. |^. 

The values of Ad? on the effective-potential loop at ^min 
and Smax are also the same because the segment of the 
curve connecting them has a positive slope. Equation 
(8c) indicates that Smax equals the lesser of | J -|- L| and 

















FIG. 4. Effective potentials ^±{S) of Eq. (14 1 for values of L, J, Si, and S 2 leading to three different sets of spin morphologies. 
The loop formed by the two curves encloses all allowed configurations for the constants listed in the legends. As in the left panel 
of Fig.|^ empty squares mark the extrema of S (Smin and 5max), empty triangles mark the extrema of ^ (^min and ^max), and 
conservation of ^ restricts the BBH spins to process along horizontal lines between the turning points S± . BBH spin precession 
can be classified into three different morphologies by the behavior of A<E> during a precession cycle: oscillation about 0 (blue 
region), circulation from —tt to n (green region), or oscillation about tt (red region). The dashed boundaries between these 
morphologies occur at values of ^ where the dotted curves cos^i = ±1 intersect the effective-potential loop, as shown by the 
empty circles. All three morphologies are present if one intersection occurs on ^+{S) and a second occurs on ^-{S) (left panel), 
oscillation of A"!* about 0 is forbidden if two intersections occur on either ^+{S) or ^-{S) (middle panel), and only oscillations 
about TT are allowed if there are no such intersections (right panel). 




FIG. 5. The (J,^) parameter space for BBHs with different minimum allowed total angular momentum Jmin. BBH spin 
morphology is shown with different colors, as indicated in the legend. The extrema ^min(>.I) and Cmax(>/) of the effective 
potentials constitute the edges of the allowed regions and are marked by solid blue (red) curves for A4> = 0 (tt). Dashed 
lines mark the boundaries between the different morphologies. The parameters q, xi, X 2 and r are chosen as in Fig. whose 
panels can be thought of as vertical (constant J) “sections” of this figure (where we suppress the S dependence). The lowest 
allowed value of ^ occurs at J = \L — Si — S 2 \ in all three panels. Three phases are present for each vertical section with 
J > \L — Si — S 2 \- This condition may either cover the entire parameter space (left panel) or leave room for additional regions 
where vertical sections include two different phases in which A4> oscillates about tt and a circulating phase in between (center 
panel) or only a single phase where the spins librate about A>1> = tt (right panel). An animated version of this figure evolving 
on the radiation-reaction time Irr is available online |55| . 
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I'S'i + S' 2 |; in the former case L and J are anti-aligned, 
while in the latter case Si and S2 are aligned. The for¬ 
mer case again requires the components of Si and S2 
perpendicular to L to be anti-aligned (A<f> = tt) but now 
the latter case requires these components to be aligned 
(A$ = 0). For values of J, L, Si, and S 2 for which this 
latter case applies, A$ will oscillate about 0 near 
and we have determined that all three morphologies are 
possible, as shown in the left panel of Fig. 

To distinguish the remaining two scenarios (whether or 
not Ad) circulates for intermediate values of ^), we must 
examine the intersections of the cos 9i = ±1 contours 
with the effective-potential loop ‘^±(5'). There can be ei¬ 
ther zero or two of such intersections. If no intersections 
occur, Ad* remains equal to tt around the entire loop just 
as it is at ^max and only oscillations about this value are 
possible, as shown in the right panel of Fig.|^ If there are 
two intersections, they must happen on the two portions 
of the loop with negative slopes (the segment connect¬ 
ing S’min and ^min and the segment connecting 5'niax and 
^max)- If both intersections happen on the same segment. 
Ad) switches from tt to 0 and back again as one traverses 
the loop from ^max to ^min resulting in the introduction 
of a circulating phase before restoring oscillations about 
TT near ^min, as seen in the middle panel of Fig. If 
the two intersections happen on different segments. Ad* 
switches to 0 at the first turning point and then to tt 
at the other leading to oscillations about 0 near ^min, as 
seen previously in the left panel of Fig. 

To summarize, the number of allowed morphologies in 
the effective-potential diagrams of Fig. [^depends on the 
magnitude of the total angular momentum J: 

1. All three phases are allowed if 

J > Si+ S 2 - L. (22) 

This condition implies S'^ax = Si + S2 and hence 
Ad>(^i„in) = 0 (Fig.| 4 j left panel). 

2. For lower values of J such that 

L-\Si-S2\<J <Si + S2-L, (23) 

Ad* will oscillate about tt near ^min and ^max and 
circulate from —tt to tt for intermediate values of ^ 
(Fig. El middle panel). The first inequality ensures 
that two (anti)aligned configurations (sin0i = 0) 
can be found, while the second prevents Ad* = 0. 

3. Finally, for 

J < min(S'i S2 — L, L — jS*! — S'2|)) ( 24 ) 

the condition sin Oi = Q cannot be satisfied and Ad* 
must oscillate about tt (Fig. |4) right panel). 

Whether these conditions can be satisfied is determined 
by the limits on J given by Eqs. In particular, Jmin = 
L — Si — S 2 is a sufficient but not necessary condition 
for all three morphologies to coexist, while J„iin = 0 is a 


necessary but not sufficient condition for the single-phase 
case. The three-phase case was considered in our Letter 
[T] and is the only allowed case at sufficiently large binary 
separations {L > Si S' 2 ). 

The J^-plane shown in Fig.[^shows all BBH spin con¬ 
figurations for fixed values of q, xi, X 2 and r at once. 
Since J and ^ are constant on the precession time tpre, the 
position of BBHs in this figure is fixed on this timescale. 
The effective-potential diagrams of Fig.E]can be thought 
of as vertical sections of Fig. at fixed J where the S di¬ 
rection has been expanded. Each panel of Fig. refers 
to a different choice of Jmin from Eqs. ( |3b[ ). Ad) can only 
oscillate about 0 if J > |L — — 5'2|. From Eq. ( [T^ , 

the limit J = |L — — S' 2 | corresponds to the lowest 

allowed value of Eor separations large enough that 
L > Si + S 2 , this configuration also corresponds to Jmin 
in which case Ad) can oscillate about 0 for all allowed 
values of J (Fig. Ej left panel). If L is sufficiently small 
to admit values of J such that J < |L — S'! — 52 !, a new 
region of the parameter space where Ad* = 0 is forbidden 
appears at small J (middle and right panels of Fig. [bj). If 
even lower values J < \Si — S 2 \ — L can be reached (i.e., 
if Jmin = 0), the leftmost region of the J^-plane does not 
even allow a circulating phase (right panel of Fig.[5]). 

The center and right panels of Fig. [5] reveal that the 
regions for which Ad) oscillates (shown in blue and red) 
are very small for L < Si + S 2 - This follows from the fact 
that these small values of the orbital angular momentum 
can only be achieved in the PN regime (r > lOM) for 
low mass ratios. Oscillation of Ad* relies upon coupling 
between the two BBH spins, and the spin S 2 becomes 
increasingly ineffective at maintaining this coupling as 
q —>■ 0 (cf. Sec. IVB below for more details). Nonethe¬ 
less, a small region of the parameter space is always oc¬ 
cupied by librating binaries as ^ approaches the resonant 
values ^min and ^max- For each value of ^ (horizontal sec¬ 
tions of Fig. [ 5 ]), one Ad) = 0 resonance and one Ad? = tt 
resonance occur at the largest (Ad* = 0) and the low¬ 
est (Ad) = tt) allowed values of J. The effective spin 
^ is therefore a good parameter to identify the resonant 
solutions, as we pointed out in Ref. m- 


D. Time dependence 


Although S fully parametrizes the precessional dy¬ 
namics, time-dependent expressions may be useful as 
well. The BBH spins obey the 2PN precession equations 

[M Eni EH 


dSi 

dt 

JS2 

dt 


1 

1 


(4 + 3 (z)L - + S2 


4 +- L- 


1 + d 


Si 


X Si, 
X S2, 


(25a) 

(25b) 


which include the quadrupole-monopole interaction com¬ 
puted in Ref. m- These equations are averaged over the 
binary’s orbital period torb and describe the evolution of 
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the spins on the precession timescale tpre- Equations ( |25[ ) 
imply that the orbit-averaged evolution of S' = |Si -I-S 2 I 
is given by: 


dt 


3(1 - g2) S 1 S 2 
2q S 

X sin 9i sin O 2 sin A$ . 


1_^\ 


(26) 


Integrating Eq. (261 yields solutions S(t), and that spec¬ 
ifies L, Si, and S 2 as functions of time through substitu¬ 
tion into Eqs. (10). Some examples of S{t) for different 
values of ^ are shown in the top panel of Fig. 

These time-dependent solutions confirm the scenario 
outlined in Sec. II B[ with S oscillating between two turn¬ 
ing points S_ and S+ at which dS/dt = 0. At these 
turning points, the three angular momenta are coplanar 
[from Eq. (26), dS/dt = 0 implies either sinA$ = 0 or 
s'mdi = 0] and the BBHs lie on the effective potentials 
{^±{S±) = ^). The spin-orbit resonances ^min and ^max 
are shown with dashed lines in Fig. [^and correspond to 
the zero-amplitude limits of the generic oscillatory so¬ 
lutions. From Eq. (26), we can define the precessional 


period r as the time needed to complete a full cycle in S, 


fS+ 

r(L,J,C)=2 / 
Js- 


dS 


\dS/dt\ ■ 


(27) 


The precession timescale tpre (27rM/77)(r/M)®/^ pro¬ 
vides an order-of-magnitude estimate for this exact pre¬ 
cessional period. The period r remains finite at the spin- 
orbit resonances ^min and ^max in much the same way 
that the period of a simple harmonic oscillator remains 
finite in the limit of small oscillations. 

The time evolution of the three angular momenta L, Si 
and S 2 is fully given by Eqs. ( [I^ and ( [^ when described 
in the non-inertial frames of Fig. However, J and L 
will generally not be confined to a plane in an inertial 
frame. The direction of J is fixed on the precession time 
scale tpre, and hence so is z. The two remaining basis 
vectors will process about the z-axis 


dx 

dt 


= fljZ X X = , 


dy 

dt 


= fl^z X y = —n^x. (28) 


The solution to these two equations gives x(t) and y{t) 
and hence L(t), Si(t), and S 2 (t) in an inertial frame from 
Eqs. (10) and ( [T^ . The orbital angular momentum L 
processes about J with frequency given by [T] 


n, = 


J 

2 

3(1 


?7" 


L2 

■d) 


2qAlAl 


1 - 


1-H 

ijM^ 

L 


2d 1 L J 

[4il-q)L^{Sl-Si) 


-{1 + q){J^ -L^ - Ar]M^L^)] 


(29) 


This equation can be derived by substituting Eqs. (25) 
and (28) into the time derivative of Eq. (§. For con¬ 
creteness, let us specify an inertial frame such that L lies 



0.0 0.5 1.0 1.5 2.0 


t/{10^M) 


FIG. 6. Time-dependent solutions for the total-spin magni¬ 
tude S (top panel) and the orbital-angular-momentum phase 
(bottom panel). We set q = 0.7, yi = O.'T y 2 = 0.9, 
r = 30M and J = 1.48M^ and integrate Eq. (26) for three 
values of ^ corresponding to the three different spin mor¬ 
phologies: A<1? oscillates about 0 (^ = 0.17, blue), circulates 
(5 = 0.25, green), and oscillates about tt (f = 0.34, red). Ini¬ 
tial conditions have been chosen such that S — S- and il?/, = 0 
at t = 0. The oscillations in S induce small wiggles in 4 ?l on 
top of a mostly linear drift. Spin-orbit resonances (horizon¬ 
tal dashed lines, top panel) correspond to configurations for 
which S is constant and can be interpreted as zero-amplitude 
limits of generic oscillatory solutions. The projections of the 
effective potentials, i.e. parametric curves [t(^)/2, 5+(5)] and 
[r(^), 5'-(^)], are shown with dotted lines. An animated ver¬ 
sion of this figure is available online . 


in the xz-plane at S' = S'-. At the point on a precession 
cycle specified by the total-spin magnitude S, the direc¬ 
tion of L is specified by the polar angles 6 l from Eq. ([^ 
and the azimuthal angle 







dS 

\dS/dt\ 


a 

2 





dS 


\dS/dt\ 


for S : S_ ^ S+ 

( 30 ) 

for S : S+ ^ S_ 
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where the two cases refer to the first and the second half 
of the precession cycle, and 


rS+ 

a{L,J,0 = 2 n. 


dS 


\dS/dt\ 


(31) 


is the total change in the azimuthal angle over a full 
precession cycle. Solutions $L(t) are shown in the bot¬ 
tom panel of Fig. The angle mainly exhibits a 
linear drift due to the leading PN order term in Eq. (251. 


Spin-spin couplings are of higher PN order and cause 
small wiggles on top of this linear drift. Binaries in spin- 
orbit resonances (■Cmin and ^max) precess at a constant 
rate flj, with all three vectors L, Si, and S 2 jointly pro¬ 
cessing about J. Just as A$ is ill defined if either of the 
Si is aligned with L (cosdi = ±1), and thus a is ill 
defined if L is aligned with J (cos^l = ±1). This occurs 
for values of J and ^ for which S- = 5'min = \ J — L\ 
or Sj^ = 5'max = J + L, corresponding to some of the 
transitions between the different classes of the evolution 
of ip' (dashed lines in Fig. [^. 

We stress here that the time-dependent expressions re¬ 
ported in this section are only valid on times t ^ r <C 
%Rj i-G. when the processional dynamics approximately 
decouples from the inspiral. This approximation breaks 
down at small separations, where the difference between 
the three timescales is smaller (cf. Sec. IIIC). 


III. PRECESSION-AVERAGED EVOLUTION 
ON THE INSPIRAL TIMESCALE 

The previous section focused on spin dynamics on the 
processional timescale. We now consider how spin pre¬ 
cession evolves as BBHs inspiral due to radiation reac¬ 
tion. Our main tool is a precession-average d equa tion to 
model the binary inspiral (derived in Sec. Ill A below) 
that will allow us to overcome numerical limitations of 
our previous analyses [SillMlIMlIinilll] and evolve BBHs 
inwards from arbitrarily large separations (Sec. IIIBl. 


This improved computational scheme relying on our new 
multi-scale analysis allows us to more efficiently “trans¬ 
fer” BBHs from the large separations where they form 
astrophysically down to the small separations relevant 
for GW detection. In Sec. |HIC we compare the results 
of our precession-averaged evolution against the standard 
integration of the merely orbit-averaged spin-precession 
equations. 


[Eq. (251 can be integrated with time steps torb <C At <C 
tpre much longer than the orbital timescale]. Radiation 
reaction can be similarly orbit averaged: 


dL 


RR 


dt 


orb 




dL 


RR 


dtp 

dt dtp / dt 


(32) 


where dL^ji/dt is the instantaneous change in the orbital 
angular momentum due to GW radiation reaction and tp 
is the true anomaly parametrizing the orbital motion. 
The flux dL^in/dt depends implicitly on both tp and the 
angular momenta L, Si, and S 2 ; the former dependence 
can be averaged over since we have analytic solutions to 
the orbital motion as function of tp, while the angular mo¬ 
menta may be held fixed, since they barely evolve over 
an orbital period. Spin precession may be calculated on 
the radiation-reaction timescale by numerically integrat¬ 
ing the coupled system of ordinary differential equations 
(ODEs) given by Eqs. (251 and (321 with the time step 
At given above. 

We derived analytic solutions to the orbit-averaged 
spin-precession equations (25) in Sec. [h] that depend on 
the magnitudes L and J that evolve on the radiation- 
reaction timescale tRR. In a similar spirit to the orbit 
averaging discussed above, we can use these solutions to 
precession average the evolution equations for L and J. 
We define the precession average of some scalar quantity 
X to be 


O /•*+ ftQ 

^ (33) 


where dS/dt is given as a function of S in Eq. (26). We 


can hold L, J, and ^ fixed on the right-hand side of this 
equation because they barely evolve during a precession 
cycle, much as we held the vectorial angular momenta 
fixed in the orbit averaging since they evolve on the longer 

timescale tpre S> torb- 

Since ^ is conserved by radiation reaction at 2.5PN or¬ 
der lilliTl, we need only find precession-averaged evolu¬ 
tion equations for L and J to evolve our spin-precession 
solutions on the radiation-reaction timescale ^rr. Since 
— L • L, dL/dt = L • dL-Q^/dt and the precession- 
averaged evolution of L is given by 


,s+ 


f) =-/ 

dt /pre ^ Js- 


dL 


RR 


dt 


dS 


orb 


\dS/dt\ ■ 


(34) 


A. Averaging the average 

In the usual PN formulation (see e.g. Ref. HU), the 
timescale hierarchy torb ^ tpre ^ ^rr is exploited to 
average the evolution equations for L, Si, and S 2 over 
the orbital period T. We already saw above how this 
orbit averaging can be used to increase the computational 
efficiency with which spin precession can be calculated 


We similarly have dJ/dt — J • dJRR/dt, but since J = 
L-l-Si -I-S 2 and GW emission does not directly affect the 
individual spins (dSi^RR/dt = 0), dJRR/dt = dL^^/dt 
and we have 


d£\ _2 / dLRR \ dS 

dt /pre Js- \ dt \dS/dt\ ■ 


(35) 


The orbit-averaged angular momentum flux 
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(dLpi,p'/(it)orb up to IPN is given by [T5] 


rfLRR \ 

dt 5 M \ r J 

2423 + 5887? M\ ^ ^ 

336 r y V / 



(36) 


Note that this expression is parallel to L and independent 


of S. Substituting this result into Eq. (35) yields 


dJ\ 

dt / pre 


dS 


ls_ ^■( d^) 


(37) 


where we used Eq. and cosOl is given in Eq. (§ as 
a function of S. Finally, Eqs. (34) and (37) together lead 
to 


dJ\ 

“ 2LJ 


(J2 + l2 _ (S'2^p^3) , (38) 


pre 


which reduces the computation of BBH spin precession 
on the radiation-reaction timescale to solving a single 
ODE [1]! Equation (|38|) is independent of the details 


of spin precession (which are encoded in (S'^)pre) and is 
also independent of the PN expansion for (dLRR/c?t)orb 
provided this is parallel to L and independent of S. As 
shown in Eq. (36), both of these conditions are satisfied 


at IPN level but break down at higher PN order. We 


address the range of validity of our approach in Sec. Ill C 


where we also perform extensive comparisons with full 
integrations of the conventional orbit-averaged equations. 

Examples of solutions to Eq. (38) are shown in Fig. 
where J is evolved from r — 10®M to r = lOM. Solutions 
J(r) are bounded at all separations by the spin-orbit res¬ 
onances ^min and ^inax whic h ext remize the magnitude 
J for each fixed ^ (cf. Sec. IIC and Fig. We per¬ 
form ODE integrations using the lsoda algorithm [62] 
as wrapped by the python module SCIPY [53]; integra¬ 
tions of Eq. (38) are numerically feasible for arbitrary 


values of<?< l,Xi ^ 1)X2 < !> and arbitrarily large 
initial separation. 

Our solutions to the spin-precession equations also de¬ 
pend on the direction J, since this defines the z-axis in 
the orthonormal frame of Fig.[2 The precession-averaged 
evolution of this direction is 



pre 


dL 


RR 


dt 


_dJ. 

orb ^ / pre 


(39) 


which is proportional to the precession average of the to¬ 
tal angular momentum radiated perpendicular to J. Al¬ 
though the vector given by the right-hand side of Eq. (|39|) 


will generally not vanish over a single precession cycle, if 
the angle a given by Eq. (311 above is not an integer mul¬ 


tiple of 27r this vector will precess about J in an inertial 
frame. This implies that J will precess in a narrow cone 


L/M^ 



r/M 


FIG. 7. Evolution of the total angular momentum magnitude 
J during the inspiral. Three binary configurations are consid¬ 
ered here: ^ = —0.5 (orange), 0 (purple) and 0.5 (green) for 
q = 0.4, xi = 0.9, X 2 = 0.8. Equation ( |38| is solved for several 
different initial conditions (solid lines, sequential colors) as the 
separation r and the angular momentum L = r]{rM^)^t‘^ de¬ 
crease. Solutions are bounded at all separations by the spin- 
orbit resonances (dotted lines) which extremize the allowed 
value of J for fixed Two of the binaries pictured here cross 
one of the resonant conditions a = 27rn (empty circles) where 
changes in the direction J are expected. The inset shows the 
same evolutions for a wider separation range. 


in an inertial frame on the radiation-reaction timescale 
remaining approximately constant [HI EH- As shown for 
some of the binaries of Fig. the condition a = 2TTn for 
integer n is indeed satisfied in generic inspirals at mean¬ 
ingful separations. Preliminary results indicate that in¬ 
teresting spin dynamics arises at these newly identified 
resonances ISD- In this paper, we restrict our attention 
to the relative orientations of the three angular momenta 
as specified by the three angles in Eqs. (20). 


B. The large-separation limit 


W e ca n gain additional physical insight by examining 
Eq. (38) in the large-separation limit L/M^ —oo. Let 
us define 


K = 


J^-L^ 

2L 


(40) 
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such that Eq. (38 1 becomes 


C. Efficient binary transfer 


dn 

dL 


{S\re 

2L2 


(41) 


The right-hand side vanishes at large separations where 
S' <C T, implying that 


= lim K 

r / M—^oo 


(42) 


is constant. This implies that n provides a more con¬ 
venient label for processing BBHs at large separations 
because it asymptotes to a constant while J diverges. At 
large separations J evolves as 


J — L) — L (2a^oo T a) , 


(43) 


as illustrated in the inset of Fig. From Eq. (40) and 
J = L -I- S one also obtains 


= lim S • L 

r / M—^OO 


(44) 


implying that k asymptotes to the projection of the total 
spin onto the orbital angular momentum. The constant 
Koo can be calculated for a binary at finite separation 
by integrating dn/dL all the way to r/M —>■ oo. This 
integration can be performed by defining u = 1 / 2 L such 
that dn/du = (S'^)pre can be integrated over a compact 
domain. 

The two constants Koo and ^ are linear combinations of 
the asymptotic values of the inner products Si • L defined 


in Eqs. (20) in the large-separation limit. The constancy 


of these inner products at large separations is also ap¬ 
parent from Eqs. (25), where the Si will precess about L 


when spin-orbit coupling dominates over spin-spin cou¬ 
pling. From Eqs. (12) and (44) one finds 


_ r f + Kao{l + q ^) 

cos 6 »ioo = hm Si • L =- „ . -r-, 

r/M—^oo ^1\Q — Q) 

(45a) 

_ 1 - of M^^-K^{l + q) 

cos 6 » 2 oo = lim Sa • L = ^ . . (45b) 

r/M^oo ^2\Q. ^ — Q) 


The terms in Eqs. (20) proportional to become in¬ 


creasingly significant at smaller separations and induce 
oscillations in 9i on the precession timescale, while the 
breakdown of the asymptotic approximation to J{L) 


given in Eq. (43) causes J (and hence 0i) to deviate on the 


radiation-reaction timescale for BBHs with different val¬ 
ues of as seen in Fig.[^ The constraints | cosdiool < 1 
and I cos daoo | < 1 define the physically allowed values of 
^ and Koo- These parameters, or equivalently dioo and 
0200 7 can be used to identify an entire BBH inspiral (as 
far as the relative orientation of the angular momenta 
is concerned) without reference to a particular separa¬ 
tion or frequency, as typically done in GW applications 


Our new precession-averaged equation for dJ/dL (38) 
can be used to efficiently “transfer” BBHs from the large 
separations at which they form astrophysically to the 
smaller separations at which the GWs they emit be¬ 
come detectable. This equation can be integrated with 
a time step ^ AT ^ much longer than the 
time step <orb At ^ tpre on which merely orbit- 
averaged equations must be integrated. This greater 
efficiency comes at the cost of no longer being able 
to keep track of the precessional phase, in much the 
same way that orbit-averaged equations do not explicitly 
evolve the orbital phase. This is not a major problem 
for population-synthesis studies however, because evolu¬ 
tion over a timescale At' will randomize the precessional 
phase, as described below. If one needs to track the pre¬ 
cessional phase below a certain separation (such as that 
corresponding to the lowest detectable GW frequencies) 
one can randomly initialize the phase at this separation 
and then employ orbit-averaged equations. The follow¬ 
ing procedure explicitly outlines how to evolve the spin 
orientations of a population of BBHs from large to small 
separations. 

1. Given a sample of BBHs specified by values of 9 , 
Xi and X 2 , choose a distribution pi(di, 02 , A$) for 
the angles that describes the spin orientations at 
an initial separation r^. This initial distribution 
is determined by the interactions between BHs and 
their astrophysical environment that lead to binary 
formation (cf. Refs. |^[57H55] on stellar-mass BHs 
and Refs. COHI!! on supermassive BBHs). 

2. Rewrite this initial distribution as a distribution 

using the relations 


S = [Sf + S 2 + 2SiS 2 (sin 01 sin 02 cos A$ 

-I- cos 6*1 008 02 )]^^^ , 

J = [L^ + + 2L{Sicos0i + S2COS02)]^^'^ , 




qSi cos 01 + S 2 cos 02 
77 M 2(1 -I- q) 


(46a) 

(46b) 

(46c) 


3. Evolve each member of the distribution fo 

a smaller separation rj using Eq. (38) for dJ/dL (^ 
remains constant). This yields a fmal distribution 

4. For each member of the distribution pf{J,^), cre¬ 
ate a distribution of values of S in the range 
S-{J,^) < S < S+{J,^) weighted by {dS/dt)~^ 
given by Eq. (26). BBHs spend less time at val¬ 
ues of S where the “velocity” dS/dt is large. This 
yields a distribution pf{S, J,^). 

5. Convert pf{S,J,^) into a distribution of final an¬ 
gles Pf {017 02 7 A<i)) using Eqs. (20) and a randomly 
chosen sign for A<i). 
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FIG. 8. Precession-averaged BBH inspirals as described in Sec. IIIC (pnrple/darker) compared to numerical integration of the 
orbit-averaged PN equations [351136j (orange/lighter). Marginalized distributions of the spin angles di, 62 , and |A4>| (rows) are 
shown at several separations along the inspirals [columns: ri = lOOOM, 500M, lOOM, SOM, and lOM]. The three initial spin 
distributions are isotropic (top panels), one aligned BH (middle panels), and Gaussian spikes (bottom panels) as described in 
Sec. |III C] The two approaches are in good agreement except for minor deviations in the distribution of A4> at r ~ lOM. We 
take q = 0.7, xi = 0-8 X 2 = 0.4 for all BBHs. An animated version of this figure is available online | 55| . 
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Examples of this binary transfer are given in Fig. [^for 
three different initial spin distributions. 


1. Isotropic sample (top panels): Both spin vectors 
are isotropically distributed (flat uncorrelated dis¬ 
tributions in cos01, cos02 and A$). 

2. One aligned BH (middle panels): One BH spin (ei¬ 
ther the spin of the primary or the spin of the sec¬ 
ondary) is aligned within 10° of the orbital angu¬ 
lar momentum, while the other spin angle 0^ has 
a flat distribution in [0°,180°]; A$ is also flat in 
[-180°, 180°]. 

3. Gaussian spikes (bottom panels): 0i and 02 have 
Gaussian distributions peaked at 45° and 135° with 
deviations of 10°; A<i) is kept flat in [—180°, 180°]. 


We evolve these distributions from = lOOOM to r/ = 
lOM and show marginalized distributions of the three 
angles 0i, 02 , and A$ at several intermediate separa¬ 
tions. An animated version of this figure can be found 
online |55j . The isotropic sample remains isotropic, as 
found previously using the orbit-averaged equations [73]. 
A greater fraction of the BBHs in the distribution with 
one aligned BH undergo a phase transition from a circu¬ 


lating to a librating morphology, as described in Sec. IV 


below and also found in previous studies with the orbit- 
averaged equations [ID]. If the angles 0^ initially have 
Gaussian distributions, these Gaussians will spread out 
as the inspiral proceeds. 

We use the BBH inspirals from = lOOOM to ry = 
lOM shown in Fig. [^ to compare the efficiency of our 
new precession-averaged approach to the integration of 
the standard - i.e., orbit-averaged - PN equations. In 
the standard approach, one must numerically integrate 
ten coupled ODEs specifying the directions of the three 
angular momenta and the magnitude of the orbital ve¬ 
locity; we use the PN equations quoted by Refs. [35ll3^ 
We implement the same 2PN spin-precession equation^ 
given by Eqs. (25) but include radiation reaction up to 
3.5PN order, as in Eq. (2.6) of Ref. [SS]. Integrations 
are performed using the same algorithm specified above 
[nH |S3] . The agreement between the two approaches is 
seen to be excellent up to r ^ 50M, and minor discrep¬ 
ancies emerge at smaller separations. 

Two approximations made in the precession-averaged 
approach may explain these discrepancies. While ^ is 
held constant throughout the inspiral in the precession- 
averaged approach (consistent with 2.5PN radiation re¬ 
action), conservation of f is not enforced in the orbit- 
averaged approach, which employs 3.5 PN radiation re¬ 
action. The largest deviations A^ in the latter approach 
are of the order 10“^°; ^ is effectively constant in the 


^ Higher-order PN corrections to the spin-precession equations 
have been computed in Refs. [ZUEZl; their implementation is 
left to future work. 


PN regime (r > lOM). Numerical-relativity simulations 
may be used to test conservation of ^ at smaller separa¬ 
tions. We have verified that additional PN corrections 


in Eq. (36), implemented in our orbit-average code up to 


3.5PN, introduce very mild corrections to the evolution 
of J: the largest variations observed in our evolutions are 
of order AJ ~ 10“^. 

The second and less reliable approximation involves 
the timescale hierarchy itself. The precession time tpre ~ 
(r/M)®/^ and radiation-reaction time tuR ^ (r/M)"^ be¬ 
come more comparable at lower separations, reducing 
the effectiveness of our quasi-adiabatic approach. The 


precession-averaging procedure defined in Eq. (33) as¬ 


sumes that quantities like L and J varying on tRR remain 
constant over a full precession cycle r, but this assump¬ 
tion will break down as the timescale hierarchy becomes 
invalid. 

Figure [^ shows that differences between the two ap¬ 
proaches are most pronounced in pj.(A$). This variable 
is the most sensitive to the precessional dynamics; pre¬ 
dictions for the angles 0i and 02 remain reasonably ac¬ 
curate even at r ~ lOM. The differences seem to aver¬ 
age out for wider distributions (top panels) but become 
more evident for more compact initial distributions (bot¬ 
tom panels). Averaging over the precessional dynamics 
prevents us from tracking the precession phase, implying 
that the two approaches will make different predictions 
for quantities (like S and A$) varying on the precession 
timescale when the initial separation is sufficiently small 
that memory of the initial phases has not been fully for¬ 
gotten. Predictions of physical quantities varying on the 
radiation-reaction timescale (like J and the precession 
morphology) will remain robust down to small separa¬ 
tions, as explored in Secs. jlV 5] and [IVC] below. 

We compare the computational efficiency of the 
precession- and orbit-averaged approaches in Fig. [^ 
Isotropic samples of 100 BBHs are transferred from large 
initial separations to a final separation Vf = lOM. 
The CPU time required by the two approaches scales dif¬ 
ferently with the initial separation. The orbit-averaged 
(OA) equations must be integrated with a time step 
shorter than the precession time, implying that the total 
number of time steps scales as 


Noa oc 


dr 


3/2 


/r/ 


^GW tpre 


(47) 


where row oc r~^ as given by the quadrupole formula 
[T9l[20]. The ratio tRR/tpre oc increases dramatically 
at large separations leading to a corresponding increase 
in the computational cost. In the precession-averaged 
(PA) approach, the integration of dJ/dL in Eq. (38) only 
requires time steps proportional to L, hence 


i-Li 

NpA oc / — log(Li) oc log(ri). 

JLr 


(48) 


The precession-averaged approach is very efficient at 
large separations because the solutions to Eq. (38) be¬ 


come very smooth in this limit as seen from Eq. (43) 
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IV. MORPHOLOGICAL PHASE TRANSITIONS 


As BBHs inspiral on the radiation-reaction timescale, 
they can transition between the spin-precession mor¬ 
phologies described in Sec. |IIC| BBH spins predom¬ 
inantly circulate at large separations but increasingly 
transition into one of the two librating morphologies 
as spin-spin coupling becomes important (Sec. 


IVAl. 


The probability of encountering one of these morpho¬ 
logical phase transitions during the inspiral depends on 
the asymmetry between the masses and the spin magni¬ 
tudes of the two BBHs (Sec. IV B[). Asymmetric binaries 


are more likely to circulate, while BBHs with comparable 
mass and spin ratios populate the librating morphologies. 
BBH spin morphologies at finite separations can be de¬ 
termined from their asymptotic spin orientations cos 0ioo 
(or equivalently ^ and Koo) as discussed in Sec. IV C 


A. Phenomenology of phase transitions 


FIG. 9. CPU time needed to evolve BBHs from an initial 
separation to a final separation r/ = lOM using our new 
precession-averaged approach (purple circles) and the stan¬ 
dard orbit-averaged approach (orange triangles). Each CPU 
time is averaged over N = 100 executions with isotropic initial 
spin orientation (flat distributions in cos0i, cos02 and Ail?). 
Dashed lines show the expected scalings: t oc for the 
orbit-averaged approach and t oc log Vi for our new precession- 
averaged approach. These computations have been performed 
on a single core of a 2013 Intel 15-3470 3.20GHz CPU. 


and Fig. Precession-averaged inspirals may even be 
computed from infinite separations through a change of 
variables to u = {2L)~^. The integrator spends most of 
the computational time at small separations, where spin 
effects - notably the numerical evaluation of S± - need 
to be tracked with high accuracy to avoid violations of 
the constraints ([^. As shown in Fig. these expected 
scalings are well reproduced by both of our codes. 


In addition to the time needed to integrate Eq. (381, 
the precession-averaged approach must generate a final 
distribution for S (step 4 above), implying that the com¬ 
putational cost does not go to zero as —)■ rj. While 
this step makes the calculation of a single BBH inspiral 
non-deterministic and more expensive, precession aver¬ 
aging effectively reduces the dimensionality of the BBH 
population during the inspiral. If the n members of this 
final distribution for S are regarded as distinct binaries, 
the total number of integrations required to produce a 
fixed number of BBHs at rj is reduced by a factor of 
n in the precession-averaged approach compared to the 
orbit-averaged approach. 


As extensively discussed in Sec. |H C| BBH spin pre¬ 
cession can be unambiguously classified into one of three 
morphologies depending on the values of q, xi, X 2 j 
r (or L), and J. While the first four of these param¬ 
eters remain constant throughout the inspiral, r and J 
evolve on the radiation-reaction timescale according to 
Eq. ( [38| ). Binaries may therefore change their preces- 
sional morphology while evolving towards merger. The 
boundaries between different morphologies (cf. Sec. H C) 
are set by the (anti)alignment condition sin0i = 0; the 
binary morphology changes whenever radiation reaction 
brings J and L to values that satisfy this condition (which 
can only occur on the effective-potential loop ^±('5'), as 
seen in Fig. |^. Figure 10 shows two examples of these 
phase transitions. At the radii rtr where phase tran¬ 
sitions occur, A$ changes discontinuously either at S- 
(left panel) or 5+ (right panel), causing the solutions 
A$(S') of Eqs. (20) to transition between the qualita¬ 
tively different shapes seen in the bottom panel of Fig. 
The BBHs in the left (right) panel evolve from the circu¬ 
lating morphology to the morphology in which Ad? oscil¬ 
lates about 0 (tt). 

A more complete phenomenology of phase transitions 
is illustrated in Fig. El The evolution of cos 9i and cos 9^ 
along the inspiral is shown for a variety of initial condi¬ 
tions cos 9iao- At each separation r, the angles 9i vary on 
the precession time within a finite range specified by the 
conditions ^ = ^±[S) (cf. Fig. |^. These envelopes vary 
on the radiation-reaction time as J evolves according to 
Eq. (38); their width shrinks to a zero as r/M —)■ oo 
according to Eqs. (45), and tends to thicken at smaller 
separations because of the in crea sing importance of terms 
proportional to in Eqs. (20). Horizontal bars above 
each panel track the binary morphologies, which we la¬ 
bel as C, LO, and Ltt for circulation, libration about 
Ad? = 0, and libration about Ad? = tt. These morpholo¬ 
gies change whenever one of the allowed ranges reach the 
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FIG. 10. Precessional solutions A<1?(S') of Eqs. (201 as J and L evolve during inspirals according to Eq. (38l. These solutions 
are colored according to the separation r/M as shown in the color bar on the right (orange/lighter for large separations 
and black/darker for small separations). Binaries in the left (right) panel transition from the circulating morphology to the 
morphology in which A<1? librates about 0 (rr) at the transition radius rtr — 152M (18.9M); separations bracketing the transition 
radius are marked with dashed lines. Parameters are set to the values indicated in the legends. An animated version of this 
figure is available online [55]. 


boundaries cos 6i = ± 1 . 

All binaries circulate at large separation because the 
angle s cos 9i and cos 6*2 are appr oximately constant 
(Sec. IIIB) and A$ from Eq. ( |20c| ) is monotonic in S', 
thus satisfying Eq. (21bI. Some binaries (leftmost panels 
of Fig. |TT| ) remain in the circulating morphology until the 
PN approximation breaks down (r = lOM). Other bi¬ 
naries undergo a single transition into a librating phase 
(middle columns of Eig. 11); A$ will oscillate about 0 
(tt) following this transition if the alignment condition 
sin^i = 0 is satisfied at S_ (S+). Since cos^i (cos 6 * 2 ) de¬ 
creases (increases) monotonically with S [cf. Eqs. (20)], 
the above conditions can be summarized as 


cos 6*1 = 1 or cos 02 = — 1 : C —y LO , (49a) 
cos 01 = — 1 or cos 02 = 1 : C — ^ Ltt . (49b) 

These phase transitions were seen in previous (orbit- 
averaged) simulations [37] and referred to as spin lock¬ 
ing, because the BBH spins locked into libration about 
the spin-orbit resonances at and ^max- As the the li¬ 
brating binaries continue to inspiral, some may transition 
back into the circulating phase, as pictured in the right¬ 
most column of Fig. [m The conditions for this second 
transition are 


cos 01 = —1 or cos 02 = 1 : LO — >• C, (50a) 
COS 01 = 1 or cos 02 = —1 : Ltt— !> C . (50b) 

As discussed further in Sec. |IVB] below, this second phase 
transition occurs in the PN regime (r > lOM) only in 
some corners of the parameter space (< 7^1 and xi ^ 
X 2 )- We have not found any additional transitions in 
the PN regime, but multiple transitions may occur at 


the smaller separations accessible to numerical-relativity 
simulations. 


B. Dependence on mass and spin asymmetry 


The asymmetry in the masses mi and spin magnitudes 
Xi determines which of the eight scenarios depicted in 
Fig. a binary will experience during its inspiral. The 
alignment conditions sin 0 i = 0 and sin 02 = 0 tend to 
be satisfied at similar values of ^ for symmetric binaries 
(g —>■ 1 and xi ~ X 2 ), shrinking the circulating (green) re¬ 
gion in the left panel of Fig. and enhancing the fraction 
of librating binaries. This point is illustrated in Fig. 12 
below, which shows the fraction of isotropic binaries in 
each of the three morphologies as functions of the binary 
separation. Each panel is computed by averaging over a 
sample of binaries isotropically distributed at large sep¬ 
arations (flat distributions in cos 0 ioo and cos 02 oo); all 
binaries in each sample share the same mass ratio and 
spin magnitudes. As the separation decreases, binaries 
transition from the circulating to librating morphologies. 
The fraction of binaries experiencing these transitions 
strongly depends on the mass ratio q. If the mass ratio 
is low (g < 0 . 6 ), most binaries remain circulating down 
to very small separations r ~ lOM. Comparable-mass 
binaries (g > 0 . 6 ) are more likely to undergo a phase 
transition in the PN regime. The typical transition ra¬ 
dius rtr at which these phase transitions occur is also very 
sensitive to the mass ratio IM1I32] ; transitions occur in 
the very late inspiral for low mass ratios while rtr can be 
as large as 10®M for g ~ 0.99. Very long evolutions are 
needed to capture all of the morphological transitions for 






































18 



FIG. 11. Evolution of the spin morphology and the allowed ranges of the spin angles 6 i over a precession cycle as functions of 
the binary separation r. Each panel shows the range of cos 6^1 (purple/darker) and cos 02 (orange/lighter) for different initial 
conditions cosOiaa- The current morphology is tracked by the horizontal bar above each panel. Morphologies are indicated as 
C (green) for circulating, LO (blue) for A<1? librating about 0, and Ltt (red) for A$ librating about tt. The morphology changes 
whenever cos0i = ±1 (vertical dashed lines). BBHs in the leftmost column do not undergo any transitions in the PN regime; 
one transition into a librating morphology occurs for BBHs in the center columns; two transitions (circulating to librating, 
librating to circulating) occur for BBHs in the rightmost column. The mass ratio and spin magnitudes are q — 0.95, Xi = 0.5, 
and X 2 = 1 in nil panels. 
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FIG. 12. The fraction / of isotropic binaries in each of the three precessional morphologies as functions of the binary separation. 
Each panel refers to different values of q, yi and X 2 as indicated in the legends. The fraction of binaries in which A>1> circulates 
(green, middle region of each panel), oscillates about 0 (blue, bottom region of each panel), or oscillates about tt (red, top 
region of each panel) is shown as the binary orbit shrinks, with dashed lines separating the different morphologies. The fraction 
of binaries in librating morphologies generally grows during the inspiral; this growth is stronger as g —>■ 1 but may stall for 
nearly equal masses and yi 7 ^X 2 , as seen in panels in the right column. 
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nearly equal-mass binaries; such long inspirals are pro¬ 
hibitively expensive in the standard orbit-averaged ap¬ 
proach (as seen in Fig. but can easily be calculated 
within our new precession-averaged formalism. 

A more extensive exploration of how BBH spin mor¬ 
phology depends on the binary parameters is shown in 
Fig. 13 and Table |Tj Isotropic distributions at r/M = oo 
are evolved down to r = lOM, where their morphologies 
are determined; as shown in the upper panel of Fig. 
these initially isotropic distributions remain isotropic at 
smaller separations. The fraction of binaries in each mor¬ 
phology at r = lOM is shown as functions of q for a grid 
of values of the spin magnitudes xi X 2 - As was al¬ 
ready seen in Fig. the likelihood of phase transitions 
depends on the mass ratio g; more librating binaries are 
found for comparable-mass BBHs at any fixed separation. 

Spin magnitudes also affect the fraction of BBHs in 
each morphology. As one moves along the diagonal of 
Fig. [^in the direction of increasing xi = A2 j a slightly 
higher fraction of binaries are found in librating mor¬ 
phologies because of increased spin-spin coupling m- 
The corner of the parameter space characterized by mass 
symmetry and spin asymmetry (g —>■ 1 and xi ^ X2) 
presents a peculiar phenomenology, as seen in the right 
panels of Fig. where the fraction of binaries in each 
morphology approaches constant values for r < lOOOM. 
This behavior can be explained by recognizing that in 
this region of parameter space binaries may undergo two 
morphological transitions in the PN regime, as seen in 
the rightmost panels of Fig. El The number of binaries 
experiencing their first phase transition from circulation 
to libration is nearly canceled by the number of binaries 
undergoing a second phase transition back to the circu¬ 
lating morphology, leading to almost constant fractions 
of binaries in each morphology. This effect also accounts 
for the kinks in the morphology fractions at g ~ 0.9 in 


the off-diagonal {xi ^ X2) panels of Fig. 13 


C. Predicting spin morphology at small separations 


We described in great detail in Sec. |IIC| how to de¬ 
termine the BBH spin morphology from the binary pa¬ 
rameters at a given separation, but astrophysical BBHs 
are often formed at much larger separations than where 
we are interested in observing them. Although BBHs 
can be efficiently evolved to smaller separations using the 
precession-averaged approach described in Sec. |HI C| we 
can in fact predict the spin morphology at a final sepa¬ 
ration Tf based solely upon the asymptotic values of 9iao 
and 02oo [or equivalently ^ and Koo according to Eqs. (45)] 
without the need to integrate dJ/dL down to ry. This 
can be achieved by recognizing that the curves in the 
cos 0ioo — cos 02 oo plane separating the final morphologies 
at Tf correspond to BBHs experiencing phase transitions 
at rf, i.e. binaries for which cos0i(r/) = ±1. These 
binaries constitute the four borders of the cos 61 — cos 62 


integrate BBHs along these borders out to r/M —)■ 00, we 
obtain four curves in the cos 0ioo — cos O 200 plane, as seen 
in Fig. These curves define regions I and II in the 
cos 0ioo — cos O 200 plane with the following boundaries: 


I. cos6»ioo=+l, 

cos0i(r/) = -1-1, 
II. cos6»ioo=-l, 

cos0i(r/) = —1, 


cos6»2oo = -1, 
cos02(r/) = -1; 
C0s6»2oo = +1, 
cos02(r/) = -1-1. 


The final morphology at r/ for each point in the 
cos 0 ioo — cos 02 oo plane is determined by whether or not 
that point is contained in the two regions: 

• Outside both region I and region II: A<I> circulates 
(no phase transitions, plain green in Fig. 14). 


• Inside region I but not region II: A$ oscillates 
about 0 (one phase transition, blue in Fig. 14). 


• Inside region II but not region I: A$ oscillates 
about TT (one phase transition, red in Fig. 14). 

• Inside both region I and region II: A$ circulates 
(two phase transitions, hatched green in Fig. 14). 


These conditions on the final morphology are consistent 
with the criteria for phase transitions given in Eqs. (49) 


and (50). Once the boundaries of regions I and II have 


plane at ry; using our expression for dJ/dL in Eq. (38) to 


been established we can determine the final morphology 
of any BBH from its initial conditions at astrophysically 
large separations without further need to integrate dJ/dL 
down to rj. A binary with spin orientations lying in the 
green, red or blue region of Fig. at large separations 
will be found with A$ circulating, oscillating about 0 or 
oscillating about tt at the end of the inspiral. 

Measuring BBH spin morphology directly offers several 
advantages over explicitly measuring the spin angles 0i, 
02 and A$. Spin morphology encodes information about 
BBH spin precession but is more robust than the spin an¬ 
gles in that it only varies on the radiation-reaction time 
(being a function of L, J, and ^). Measurement of only 
the two angles 0 i and 02 at small separations constrains 
neither the morphology at small separations nor the ini¬ 
tial conditions at large separations, as can be seen from 
the scatter points in Fig. which show an isotropic 
sample of binaries at rf. Points corresponding to the cir¬ 
culating and both librating morphologies lie right on top 
of each other in this plot, evidence of both the impor¬ 
tance of the third angle A<I> and the large oscillations in 
0i at small separations seen in Fig. El By contrast, spin 
morphology is a direct memory of a BBH’s initial posi¬ 
tion in the cos 9iao — cos 6*200 plane, as seen in Fig. 
Astrophysical scenarios of BBH formation can favor some 
regions in this plane over others [40j, implying that GW 
observations of spin morphology can constrain BBH for¬ 
mation im. 
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FIG. 13. The fraction / of isotropic BBHs for which A<1> circulates (green, middle region), oscillates about 0 (blue, bottom 
region), or oscillates about tt (red, top region) at a binary separation r = lOM as functions of the mass ratio q. Dashed lines 
separate the different morphologies. Each panel corresponds to a different value of Xi (columns) and X 2 (rows). The fraction 
of BBHs in librating morphologies increases as the mass asymmetry decreases (? —^ 1). For nearly equal masses {q > 0.9), 
asymmetry in the spin magnitudes increases the fraction of binaries in the circulating morphology as can be seen by comparing 
panels on and off of the diagonal. Some data used in this plot are listed in Table [I] The website contains an animated 
version of this figure, where the panels are shown at decreasing binary separations. 
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FIG. 14. Spin morphologies at r/ = lOM as functions of the asymptotic values of the spin angles 9ioo- The mass ratio q and 
spin magnitudes Xi for each panel are indicated in the legends. Evolving BBHs along the four lines cos^^ = ±1 at r/ out to 
r/M —>■ 00 using our new precession-averaged approach yields the dashed curves separating the different final morphologies: A$ 
oscillates about 0 (blue), oscillates about rr (red), circulates without ever having experienced a phase transition (plain green), 
or circulates after having experienced a phase transition to libration and then a second phase transition back to circulation 
(hatched green). The morphology within each region defined by the dashed boundaries is determined by which of the conditions 
cosOi = ±1 these boundaries satisfy, as described in Sec. |IV C] The points show the locations of binaries in the cos0i — cos 6^2 
plane at r/ and are colored by their morphology at that separation [A4> oscillates about 0 (blue circles), oscillates about tt 
( green squares), or circulates (red triangles)]. Because morphology depends on A4> in addition to 9i and 02 at finite separation, 
the projection onto the cos 9i — cos 02 plane can lead points of different morphologies to occur at the same positions, particularly 
for comparable-mass binaries g ~ 1 where the 0i’s oscillate with greater amplitude. The website contains an animated 
version of this figure in which r/ evolves. 


V. DISCUSSION 


BBHs evolve on three distinct timescales: the orbital 
time torb, the precession time tpre, and the radiation- 
reaction time tRR. In the PN regime (r >> Vg), these 
timescales obey a strict hierarchy: torb ^pre ^ ^rr- 
All of the parameters needed to describe BBHs evolve 
on a distinct timescale: the vectorial binary separation 
r on torb, the angular-momentum directions L and 
on tpre, and the orbital-angular-momentum magnitude L 
and total angular momentum J on Irr. The mass ratio q 
and spin magnitudes Si remain constant throughout the 
inspiral. Expanding on our previous Letter [1], we exploit 
this timescale hierarchy and conservation of the projected 
effective spin ^ [Ml [57] throughout the inspiral to solve 


the orbit-averaged 2PN equations of BBH spin precession 
given by Eq. (25). The solutions given by Eq. (20) for 
the three angles 0i, 02, and A$ that specify the relative 
orientations of L, Si, and S2 are remarkably simple and 
are given parametrically in terms of a single variable, the 
total-spin magnitude S', that evolves on tpre. 


These solutions fully determine how the relative orien¬ 
tations of the three angular momenta evolve over a pre¬ 
cession cycle as S oscillates back and forth between ex¬ 
trema S±. We find that spin precession can be classified 
into three distinct morphologies depending on whether 
A<1> oscillates about 0, oscillates about tt, or circulates 
through the full range [— tt, -|-7r] over a precession cycle. 
For BBHs with a given mass ratio q and spin magni¬ 
tudes Si, the precessional morphology at a binary sep- 
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Xi = 0.2 xi = 0.4 xi = 0.6 Xi = 0.8 Xi = 1 
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TABLE I . Fractions of isotropic BBHs in each of the three precessional morphologies ( LO : A4> oscillates about 0 , C : A$ 
circulates, Ltt : A<f> oscillates about tt ) at r = lOM as shown in Fig. |13[ For a grid of values in Xi (columns), X 2 (rows) and, 
q (first column in each mini-table), we report the fraction of binaries in each morphology. The sum of the three fractions may 
differ from unity because of rounding errors. 
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aration r is determined by J and implying that the 
morphology only evolves on the radiation-reaction time 
tRR- Spin-orbit coupling dominates over the higher-PN- 
order spin-spin coupling at large separations implying 
that all BBHs formed at such large separations begin 
in the circulating morphology. Since ^ is constant to 
high accuracy throughout the inspiral, evolving our so¬ 
lutions (20) and their associated morphology to smaller 
separations (lower values of L) only requires an expres¬ 
sion for dJ/dL due to radiation reaction. All previous 
studies of radiation reaction have relied on orbit-averaged 
expressions for dh^iYi/dt that must be integrated numer¬ 
ically with time steps At < tp^e- Our new solutions (20) 
allow us to precession average these expressions to de¬ 
rive Eq. (38) for dJ/dL that can be integrated with a 
time step tpre At' < tRR. The computational cost of 
calculating inspirals from an initial separation in our 
new precession-averaged approach scales as logr^, lead¬ 
ing to vast savings over the traditionally orbit-averaged 
approach (which scales as for the large initial sep¬ 

arations relevant to astrophysical BBH formation. 


Using our new expression for dJ/dL^ we can evolve 
our initially circulating BBHs to smaller separations, 
where they may experience a phase transition to one of 
the two librating morphologies. Some of these librating 
BBHs may subsequently undergo a second phase transi¬ 
tion back to circulation before reaching a binary separa¬ 
tion r = lOM below which the PN approximation itself 
begins to break down. Our precession-averaged calcu¬ 
lation of the inspiral agrees well with the orbit-averaged 
approach down to nearly this separation where small dis¬ 
crepancies appear because of dynamically generated in¬ 
homogeneity in the precessional phase as the timescale 
hierarchy fails. Unlike the angles 0i, 02 and A$, that 
vary rapidly on the precession time at small separations, 
the precession morphology at small separations is directly 
determined by the asymptotic values Qiao of these angles 
at large separations, providing a memory of BBH forma¬ 


tion potentially accessible to GW detectors. 

Although this work focuses on BBH spin precession, 
our analysis also facilitates the calculation and inter¬ 
pretation of GW signals. Fast templates suitable for 
GW detection and parameter estimation are being devel¬ 
oped using our new precessional solutions and precession- 
averaged equation for radiation reaction |78j . The in¬ 
sights underpinning our approach (most notably the use 
of a hierarchical coordinate system that better respects 
the separation of timescales intrinsic in the binary dy¬ 
namics) are also helping us to assess whether the preces¬ 
sional morphology of BBHs in spin-orbit resonances can 
be reliably identified in the context of full GW param¬ 
eter estimation m- Preliminary results indicate that 
BBH spin orientations can be significantly constrained 
at realistic signal-to-noise ratios, suggesting that obser¬ 
vations of BBH spin precession as described in this work 
may soon provide a new window into the astrophysical 
origins of BBHs and general relativity itself. 
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